library(dplyr)
library(tidyr)
set.seed(12345)

if (!file.exists("output")) {
  dir.create("output")
}

# load data
source("HNC_metadata_tidy.R")
ancestry <- read.csv("Supplementary_Table_5_data_ancestry.csv") %>% select(-country)
# CN clusters generated in Figure_7c.R
CN_clusters <- read.csv("output/CN_clusters.csv")

data <- data %>%
  left_join(ancestry, by = "donor_id") %>%
  left_join(CN_clusters, by = "donor_id") %>%
  filter(!is.na(HNC_clusters))

var1 <- "HNC_clusters"
var2 <- c("tobacco_ever", "alcohol_ever", "tob_alc", "hpv_pos_opc", "bmi_cat", "oral_cat", "hotdrinks", "subsite", "sex", "age_group", "stage", "region", "country", "predicted_ancestry")

fisher_table <- NULL
for (f in var1) {
  for (v in var2) {
    if (exists("test") == T) {
      rm(test)
    }

    samplesize <- data %>%
      group_by(data[, f], data[, v]) %>%
      tally() %>%
      mutate(
        perc = prop.table(n) * 100,
        perc = paste0(n, " (", round(perc, 1), ")"), .keep = "unused"
      ) %>%
      mutate(variable2 = v, .before = everything()) %>%
      spread(`data[, f]`, perc, fill = 0) %>%
      rename(levels = `data[, v]`) %>%
      filter(!is.na(levels))

    t <- table(data[, v], data[, f])
    if (ncol(t) < 2 | nrow(t) < 2) {
      print(paste("Test for v1", f, "v2", v, "has been removed due to <2 categories"))
    } else {
      test <- tryCatch(
        {
          fisher.test(t)
        },
        error = function(e) {
          NULL
        }
      )

      if (is.null(test)) {
        print(paste0("Performing Fisher test with simulated p value for v1 ", f, " v2 ", v))
        test <- fisher.test(t, simulate.p.value = T)
      }

      result <- data.frame(
        variable1 = f,
        variable2 = v,
        fisher_p = test$p.value
      ) %>%
        right_join(samplesize, by = "variable2", relationship = "many-to-many")
      fisher_table <- rbind(fisher_table, result)
    }
  }
}

fisher_sig <- fisher_table %>%
  select(variable1, variable2, fisher_p) %>%
  unique() %>%
  group_by(variable1) %>%
  mutate(
    variable2 = factor(variable2, levels = variable2),
    q_val = p.adjust(fisher_p, "fdr")
  ) %>%
  ungroup() %>%
  select(variable2, q_val)

fisher_table <- fisher_table %>%
  left_join(fisher_sig, by = "variable2") %>%
  select(-fisher_p, -q_val, everything(), fisher_p, q_val)

write.csv(fisher_table, "output/Supplementary_Table_17.csv", row.names = F)
